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1 Introduction 

Particle multiplicity distribution is one of the first measurements to be taken at the 
CERN Large Hadron Collider (LHC). The upcoming data on its energy, centrality, and 
rapidity dependence are expected to discriminate various (integrated and unintegrated) 
parton distributions, which are basic qualities analyzing high-energy reactions on the 
parton level. The Bjorken variable x of gluons may reach very small values at LHC 
energies. Therefore, the nonlinear corrections of the initial gluon correlations to the QCD 
evolution equations should be considered in any available parton distributions of the LHC 
physics [1]. 

One of the striking predictions of nonlinear QCD evolution equations is the saturation 
solution of the Jalilian-Marian-Iancu-McLerran- Weigert-Leonidov-Kovner ( JIMWLK) equa- 
tion [2], where the unintegrated gluon distribution is absolutely flat in kt-space at kt < Qs', 
Qs is the saturation scale and the relating phenomenon is called the color glass conden- 
sate (CGC). The above saturation solution was proved in the numerical solutions of the 
JIMWLK equation [3]. Instead of the complicated JIMWLK equation, which is equiv- 
alent to an infinite, coupled hierarchy of evolution equations (the Balitsky hierarchies), 
much of the physics content has been studied with the help of approximations, such as the 
Balitsky-Kovchegov (BK) equation [4], or outright models that capture the main features 
as they affect specific phenomena, such as the Kharzeev-Levin (KL) model [5] does for 
the present context. 

The KL model mimics a possible saturation solution of the JIMWK equation, and it as- 
sumes that the scale Q'^{x) = Qq{xq/x)^. Under some additional assumptions, Kharzeev, 
Levin, and Nardi (KLN) use this model successfully to explain some of the climactical data 



probed at the BNL Relativistic Heavy-Ion Collider (RHIC) and to predict LHC physics 
[6]. The BK equation is regarded as the leading order approximation of the JIMWLK 
equation and preserves the saturation and scaling features. Albacete uses the BK equation 
to predict the pseudorapidity density of charged particles produced in central Pb — Pb 
collisions at LHC [7]. 

However, both the KL model and BK equation work at a very small x range, (for say, 
X < Xq, Xq ^ 10~^), where the saturation begins work. The x-dependence of the gluon 
density at a; > Xq in these works [6,7] are simply fixed to be (1 — x)'^, which is un-evolution. 
As we have pointed out, the parton distributions at intermediate and lager x influence 
the shape of the rapidity spectrum [8]. In fact, Szczurek uses the same KL model but 
with different factors ~ (1 — x)^^^ and obtains a narrower pseudo-rapidity distribution 
[9]. Anyway, a fixed form ~ (1 — x)'^ for an available gluon distribution is too rough. 
Obviously, as a complement to the above-mentioned saturation models, it is necessary to 
consider the parton distributions, which are well defined in a broad kinematic range. 

Besides, some works [6,7] have used a single production mechanism, i.e., gluon-gluon 
fusion gg — )■ g, which is proposed by Gribov, Levin and Ryskin (GLR) [10] to predict 
particle production. However, the single-particle inclusive spectrum shows that its rapid- 
ity distribution has three distinct regions: a central region with two (project and target) 
fragmentation regions [11]. The GLR model is expected to dominate the processes at the 
central region. Therefore, other collision dynamics in the fragmentation regions should be 
considered. There are several two-component models for hadron collisions. For example, 
the HIJING model [12] is such a two-component model. This model uses the parameter- 
ized integrated parton distributions to compute multiple minijet production and incor- 



porates the Lund string model [13] to model soft beam jet fragmentation. However, the 
production of minijets at the central region should be described by the unintegrated gluon 
distribution rather than the integrated gluon density [9]. Moreover, the string model is 
irrelevant to the parton distributions. In the same model, the hard and soft components 
based on different basic physical parameters may lose some interesting information about 
heavy-ion collisions. 

In this work, we use a two-component model and a set of consistent parton distri- 
butions to improve the above-mentioned situation of the present models. Basically, we 
use two complementary production mechanisms: the hard gluon-gluon fusion [10] and 
the soft quark recombination [14]. Our picture is as follows. At sufficiently high-energy 
hadron-hadron collisions, particles produced in hard gluon-gluon fusions are distributed 
in a region around midrapidity, and these initial gluons are described by the unintegrated 
gluon distribution in both colliding hadrons. On the other hand, the valence quarks tend 
to fly through with their original integrated distributions and they are hadronized by 
recombining with additional low- a; sea quarks from the central region. The resulting soft 
recombined particles dominate the fragmentation region. We shall present the related 
formula of two production mechanisms in Sec. II. Because gluon fusion and quark recom- 
bination use unintegrated and integrated parton distributions, respectively, a set of these 
parton distributions, which are defined in a broad kinematic range are necessary. In our 
previous work [8] , we proposed such distributions in protons and heavy nuclei by using a 
modified Dokshitzer-Gribov- Lipatov-Altarelli-Parisi (MD-DGLAP) equation [15], which 
incorporates the shadowing and antishadowing corrections. Unlike the JIMWLK and 
BK equations, MD-DGLAP equation works in a broad pre-saturation range. We shall 



use these parton distributions to explain the particle multiplicity distributions in UA5 
and RHIC data, then we will predict the particle multiplicity distributions at LHC in 
Sees. Ill and IV. Unfortunately, we can not theoretically predict the energy dependence 
of the normalization of the hadronic distributions, which is fixed using the entirely phe- 
nomenological parametrization of existing data in this work. This simple method may 
need modification after we obtain the LHC data at midrapidity. 

In Sec. V we discuss some interesting properties of rapidity distributions. Limiting 
fragmentation [16] and rapidity plateau [17] are characteristics of two components of par- 
ticle production. We find that the limiting fragmentation hypothesis, which generally 
appears in present data of hadron collisions is partly violated if the observations are over 
a wide range between the RHIC and LHC energies. An explanation of limiting frag- 
mentation and its violation in partonic picture is given. On the other hand, we propose 
that a possible quark gluon plasma (QGP) effect may deform the shape of the central 
plateau. We also discuss the nuclear shadowing effects in heavy-ion collisions. The com- 
parisons of our predictions with those of saturation models are given in this section. The 
above-mentioned properties of particle multiplicity distributions in the partonic models 
of hadron collisions, except for their height at midrapidity, are dominated by the parton 
distributions. Therefore, the results potentially tell us how parton distributions evolve at 
high energies. The last section is a short summary of this study. 



2 Two-component model 

The single-particle inclusive spectrum shows a rapidity distribution with three distinct 
regions: a central region and two (project and target) fragmentation regions. We assume 
that the hadrons produced in the central region (small x and large kt) are produced 
from the hadronization of the gluons in the gg ^ g mechanism [10], while the particles 
in the fragmentation region are formed by the valence quarks according to the quark 
recombination model [14]. The related formulas are summarized as follows. 

Component I. The cross section for inclusive gluon production in pp ^ g through the 
gluonic mechanism gg ^ g ai sufficiently high energy reads [10] 
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where fi = max{klt, k^t^PtJ , kft = l{pt,g + qt,gf and kit = \{pt,g - gt,^)^ The rapidity 
y of the produced gluon in the center-of-mass frame of p — p collisions is defined by 



Xi/2 = ^ ■ exp{±y); (2) 



F^{x,k1,p'l) is the two scale unintegrated gluon distribution in the proton. A general 
relation between integrated and unintegrated parton distributions is 

/ dk^F^ix, kl, /i^) = xap{x, /i^), (3) 
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where ap{x,iJ?) = Vp{x,fi'^), Sp^x^jj?) and gp{x,fj.'^) imply the integrated valence quark, 
sea quark and gluon distributions in the proton. 



In experiments a good identification of particles is to measure pseudorapidity. Tlie 
relation between the rapidity y and pseudorapidity t] for massive particles is 
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where ^T^e// is the typical invariant mass of the gluon mini-jet. 

For avoiding the complicate hadronization dynamics, similar to Ref. [9], we use local 
parton-hadron duahty, i.e., the rapidity distribution of particles is identical to the rapidity 
distribution of gluons: rjg = rjh = t]. Thus, the pseudorapidity density of produced charged 
particles in p — p collisions is given by 
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where the Jacobian is 
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and we neglect the fragmentation function D{z,pt^h)- Corrections to the kinematics due 
to the hadron mass are considered by replacing pi -^ pi + 'miff in the evaluation of a;i/2- 
Assuming pions in p~p collisions are produced via p-resonance, we take rrieff = 770MeV. 
Component II. According to the quark recombination model [14], the valence quarks 
of incident proton tend to fly through the central region with their original momentum 



fraction. These valence quarks recombine with lower pt antiquarks and produce the out- 
going hadrons in the fragmentation region. The quark recombination model has explained 
successfully the meson inclusive distributions and the leading particle effects in the frag- 
mentation region below RHIC energies [18]. 

The cross section for inclusive pion production in p — p collisions in the quark recom- 
bination model is 
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where dspi^x^pj) is the distribution of additional sea quarks in the central region and we 
assume that it has the form like Sp{x,pf). 

We introduce the rapidity for pions in the recombination processes as 
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Thus, we have 
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and they lead to 
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where we integrate over transverse momenta in Eq. (11) at pt < IGeV, since the recom- 
bination model works at lower pt range. 

Summing the contributions of two components, we have the total distribution 

dN^^dl^^dNlL^^ 
dr] drj drj 

The production mechanisms (I) and (II) use unintegrated and integrated parton dis- 
tributions, respectively. In particularly, the gluon momentum fraction in Eqs. (1) and 
(7) contain both smaller and larger x regions. Therefore, a set of consistent integrated 
and unintegrated parton distributions in protons and heavy nuclei, which are defined 
in a broad kinematic region, are necessary. Fortunately, such parton distributions are 
proposed in Ref. [8], where the integrated parton distributions are evolved by using a 
modified DGLAP equation [15] in a whole pre-saturation region; while the unintegrated 
parton distributions are obtained directly from these integrated parton distributions using 
the Kimber, Martin and Ryskin (KMR) scheme [19]. We shall use these parton distribu- 
tions to predict the particle multiplicity distributions. 



3 Proton-proton collisions 

We calculate the pion distributions in p — p collisions using the two component model 
and compare the results with the data. We need the values of the total inelastic cross 
section o"i„, which is included by the coefficient c{^/s) of Eq. (5). This is rather compli- 
cated at the parton level since it contains the nonperturbative information. In this work, 
we use the midrapidity density for p — p{p) collisions to estimate the values of c{^/s). The 
former has been parameterized by the UA5 [20] and CDF [21] collaborations as 



dNp_p 



= 2.5 -0.25 In s + 0.023 In^s. (13) 
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This purely empirical parametrization is fitted in a broad range from ^/s = 15GeV to 

l.STeV and we extend it to the LHC energies. The resulting c{^/s) is shown in Fig.l. 

Figure 2 shows our results for dNp^p/drj bellow LHC energies with rrieff = 770MeV 
(solid curves). The data are taken from Ref. [20]. For comparison, we also draw the 
distributions with me// = 0. We find that the shape of the central rapidity plateau 
relates sensitively to the value of parameter me//. As shown in Fig. 2, corresponding to 
me// = 770MeV, the full rapidity plateau has two peaks. During the reduction of the 
nieff value, the rapidity plateau is fiattened at \ri\ < 2 and even disappears. 

To illustrate the contributions from two production components in our model, in Figs. 3 
and 4 we use dashed and dotted curves to indicate the contributions from the gluon fusion 
and quark recombination models, respectively. Although the saturation models [6,7] have 
used a single gluon fusion mechanism to reproduce these data, they both take a fixed 
gluon distribution ~ (1 — a;)'^ at the pre-saturation range. Now we use a reasonable gluon 
distribution instead of (1 — a;)^. One can find that the resulting dashed curves are narrower 
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than the sohd curves (which are consistent with the data in Fig. 2) in Figs. 3 and 4, and 
this imphes that an additional contribution from the fragmentation regions is necessary. 

We predict the pion distributions in p — p colhsions at LHC energies in Figs. 5-9. Al- 
though the contributions of the quark recombination model (dotted curves) are generally 
smaller than that of the gluon-gluon fusion model (dashed curves), the contributions of 
the quark recombination still can not be neglected in the fragmentation region. 

Several partonic models have predicted the pion distributions in p — p collisions at 
LHC energies. It is interest to compare our results with them. Figure 10 presents the pre- 
dictions of KLN work [6] (dashed curve) and the comparison with our result (solid curve). 
Figure 11 compares the comparisons of our prediction with those of the PYTHIA model, 
which is based on string fragmentation mechanism [22], and PHOJET model, which uses 
a Pomeron exchange [23]. Figure 12 is a similar comparison with the ultrarelativistic 
quantum molecular dynamics (UrQMD) model [24]. We think that the comparisons with 
different models can provide useful knowledge about the unintegrated gluon distribution 
in the proton and a correct picture of hadron collisions. 
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4 Nucleus-nucleus collisions 

The nucleus-nucleus collisions are much more complicated than p — p collisions. The 
high multiplicities in heavy-ion collisions typically arise from the large number of nucleon- 
nucleon collisions. In the analysis of heavy-ion collision data at highly relativistic energies, 
two parameters which characterize the influence of nuclear geometry are used [25]: (1) the 
number of participating nucleons Npart, which depends on the collision geometry, and (2) 
the number of binary nucleon-nucleon collisions Ncoii, or the average struck number of each 
participating nucleon as it passes through the oncoming nucleus, i7 = Ncoa/{0.5Npart)- In 
nuclear collisions, the soft (or hard) component is proportional to the number of partici- 
pants Npart (or the number of binary colhsions Ncou) [26]. Essentially, we write 

< Npart > dr] y * drjdp'i Oin J * drjdp'^ 

where < ... > is an average value in a giving central cut and we only consider — 6% cut in 
this work; da a- a implies that the parton distributions are the nuclear parton distributions. 
According to the geometric approach in Ref.[27], we take the mean number of participants 
Npart = 339 in Au — Au collisions at ^/s = 130 and 200GeV, and Npart = 369 in Pb — Pb 
collisions at y/s = 5.5TeV. 

The value of u contain the knowledge about the interaction between two collided 
nuclei. Glauber modeling in high-energy nuclear collisions [28] has pointed out that the 
number of collisions roughly is 



< N,M >oc< N^i^^t >=\< N^lt > 



(15) 
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On the other hand, the RHIC data [29] present a slow increase of the coefficient A with 
energy ^/s. Using these data about dNAu-Au/dri\ri=o at ^/s = 130 — 200GeV, we take a 
best ffiting: A = 29 In ^/s. 

Now we can "predict" the whole distributions in Au — Au collisions at ^/s = 130 and 
200GeV. The results are shown in Figs. 13 and 14. Where we keep temporarily the value 
of parameter nieff = 770MeV. The dashed and dotted curves are the contributions of the 
gluon-gluon fusion model and quark recombination model, respectively. There is a small 
deviation from the data at \ri\ > 5 since the Fermi motion contributions [30] are neglected 
in our considerations. 

Considering the experimental errors, the inconsistency between the theoretical curves 
and the RHIC data is still visible, and it suggests that some factors failed in our above 
considerations. In Figs. 15 and 16, we reduce the value of Tn^ff. A best fit (solid 
curves) requests the parameters me// = 400MeV. Compared with our previous result of 
me// = 770MeV for p — p collisions in Sec. Ill, we find that the reduction of Trieff at 
RHIC Au — Au collisions is possible. 

Finally, we calculate the pseudorapidity distributions in Pb — Ph collisions at ^/s = 
5.5TeV. The differences between lead and gold are neglected. Our results are shown by 
sohd curve in Fig. 17, where we take me// = 400GeV; the dashed, and dotted curves corre- 
spond to the contributions from gluon-gluon fusion and quark recombination, respectively. 
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5 Discussions 

(i) Limiting fragmentation, is it universal or not? 

To separate the trivial kinematic broadening of the dN/drj distribution from more 
interesting dynamics, the collision data at different energies are viewed in the rest frame 
of one of the colliding target. Such distributions lead to a striking universality of multi- 
particle production-limiting fragmentation. This hypothesis states that, at high enough 
collision energy, when effectively viewed in the target rest frame, dN/drj' exhibits lon- 
gitudinal scaling and becomes independent of energy in a region around rj' ~ 0, where 
f]' = f] — ybeam, [ybeam = ^^{y/s / m n)] . The hypothcsis of limiting fragmentation in high 
energy hadron-hadron collisions was first suggested in Ref. [16]. From a phenomenological 
view, the projectile hadron, when seen in the frame of the target, is Lorentz-contracted 
into a very narrow strongly-interacting pancake which passes through the target, assum- 
ing that the total hadronic cross sections would become constant at large center-of-mass 
energy. If this occurred, the excitation and break-up of a hadron would be independent of 
the center-of-mass energy and distributions in the fragmentation region would approach 
a limiting curve. We know that the total hadronic cross-sections are not constant at 
high energies, therefore, limiting fragmentation should fail. However, limiting fragmenta- 
tion has been observed in a wider region, even extending nearly to midrapidity, and it is 
referred to as extended longitudinal scaling [25]. 

From the partonic point of view, longitudinal scaling in hadron collisions relates to 
Bjorken scaling of the parton distributions and the production dynamics. An interesting 
question is whether two component models can keep the limiting fragmentation curve. To 
answer this question, we plot the shifted pseudorapidity distributions in central Au — Au 
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collisions at ^/s = ISOGeV and 200GeV in Fig. 18. The distributions are scaled by Npart/'^ 
to remove the effect of the different number of nucleons participating in the collisions. We 
find longitudinal scaling (energy independence) over more than three units of rapidity, 
extending nearly to midrapidity and it is consistent with the RHIC data [31]. In Fig. 
19, we present the contributions only from the gluon fusion mechanism, where limiting 
fragmentation still holds at ?]' > — 1. 

However, we compare the similar distributions including Ph — Pb collisions at ^/s = 
5.5TeV (dotted curve)in Fig. 20. We find a smaller deviation from the limiting frag- 
mentation limit at 1]' < 0, although the distributions at r^' > still keep longitudinal 
scaling, where it is dominated by the quark recombination mechanism. Figure 21 shows 
the comparisons of the contributions from the gluon-gluon fusions. The results indicate 
the deviation from limiting fragmentation origins from the gluon-gluon fusion mechanism. 

Back to the p — p collisions. In Fig. 22 we plot our predicted curves in p — p collisions 
from ^/s = 130GeV to 14TeV with t]'. We find that a similar deviation from limiting 
fragmentation exists at —2 < r^' < if the energies across over a big range between RHIC 
and LHC. 

We noted that a different deviation from limiting fragmentation at the LHC en- 
ergy is also predicted by using another kind of parton distributions, i.e., the McLerran- 
Venugopalan distributions in Ref. [32], where the same gluon-gluon fusion mechanism is 
used. Therefore, investigation of a possible deviation from the limiting curve will provide 
insight into the evolution equations for high-energy QCD, although the possible larger 
systematic errors in experiments may hide the deviation. 

We try to understand limiting fragmentation and its violation from the partonic pic- 
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ture. The gluon distributions in Eq.(l) are really irrelevant to the interaction energy y/s 
in the parton model [17]. A possible relation of Eq. (1) with the interaction energy is that 
the kinematic ranges of Bjorken variables a;i/2 are y^-dependent. To illustrate that, we 
draw the kinematic ranges of two multiplying distributions at three different energies us- 
ing Eq.(2) in Fig. 23. For example, we fix pt = O.SGeV" in Eq.(2) and take mf,ff = 0, thus 
we have 1] = y. We can find that at y = y^m (or y' = y — ybim = 0), an extremely small 

Xl = Xi^small (or X2 = X2,small) alwayS combines with a larger X2 = X2,large (or Xi^iarge)- 

Besides, 



rriN 



is independent of ^/s, and 



Xl^small — X2,small — • U'J 



Thus, we have 



Fg[Xijarge, Qt, Pt)Fg{x2,small, Qt, Pt)\y'=0,y/^=200GeV 
— Fg{xi^large,qt,Pt)Fg{x2,smaU,qt,Pt)\y'=0,^=l30GeVi (18) 

since X2,smaii(}V = ISOGeV) ~ X2,smaii{y/s = 200G'eV"). On the other hand, the difference 
between the parameters c{y/s) at ^i = 130 and 200GeV is small. Consequently, we have 
limiting fragmentation near y' = in the gluon-gluon fusion processes in Fig. 20. 

However, an obvious difference exists near y' = ii we compare the results at the 
RHIC and LHC energies, i.e.. 
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Fg{Xi^large, Qt, Pt)Fg{x2,small, Qt, Pt)\y'=0,^=200GeV 
^ Fg{xi^large, qt,Pt)Fg{x2,smalh qt,Pt)\y'=0,y^=5.5TeVi (19) 

since X2,smaui\^ = 5.5TeV) « X2,smaiii\/s = 200GeV). Although the decreasing c(^/s) 
with increasing ^/s almost compensates for the difference in Eq.(19), o"j„ ~ y/s and 
Fg ~ ^/s belong to really different dynamics; therefore, the results in Fig. 21 show a 
deviation from limiting fragmentation. 

It is different from the gluon-gluon fusion mechanism; the quark recombination model 
naturally satisfies limiting fragmentation. The reason is as follows. For a given y' = 



Inx + \n{mj\f/Jml. +Pt), the kinematic ranges of Eq. (7) are Xi E [x,0] and X2 G [0,x], 
which are irrelevant to the energy ^/s. Therefore, the resulting rapidity distributions have 
limiting fragmentation, as we have shown in Fig. 20. 

(ii) The possible deformation of central rapidity plateau. 

A general picture of two components of particle production in hadron-hadron collisions 
predicts two types of ranges for the distributions of final-state particles: except for limiting 
fragmentation at the fragmentation ranges, particles near midrapidity in the center-of- 
mass frame were expected to form a rapidity plateau [17]. A narrow plateau appears in 
our distributions, its height (but not the width) grows with energy. We point out that 
the structure of the central plateau relates to the value of parameter rrieff rather than 
the parton distributions: a larger value of me// may structure a plateau with double 
peaks, while it is flattened with decreasing Trieff, and the plateau even disappears when 
me// -^ 0. 

Comparing a best fitting me// = 770 MeV in p — p collisions at ^/s = 200GeV (see Fig. 
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2) with that of Trieff = 400MeV in Au — Au colhsions at the same energy (see Fig. 16), 
the reduction of Trieff is possible in ultra relativistic heavy-ion collisions. An interesting 
question is: whether this is a QGP effect in the RHIC data. 

One of the most exciting research areas in RHIC collisions is to find a new matter, 
i.e., the QGP, in which quarks and gluons are no longer confined to volumes of hadronic 
dimensions and hadron masses are reduced under the QCD phase transition. One of the 
conditions in forming the QGP is that a sizable fraction of the initial kinetic energy creates 
many thousands of particles in a limited volume. Therefore, we consider that the QGP 
is formed after gluon-gluon collision in the central region in ultrarelativistic heavy-ion 
collisions. The effective mass meff in Eq. (4) will be lowered if part of the minijets go 
through the QGP region, since that is where chiral symmetry is restored. We regard this 
decreasing value of meff as a QGP effect partly due to the restoration of chiral symmetry. 
Of course, since large errors occur with the RHIC data in the region rj < 3 [31], more 
accurate data are necessary. However, if the above-mentioned QGP corrections to the 
RHIC data are true, one can expect a more obvious effect in Pb — Pb collisions at LHC. 
In Fig.24, we plot such effects in the central plateau of Pb — Pb collisions at ^/s = 5.5 
TeV where the values rrieff are reduced from 400 MeV to 0. 

(iii) How large are the nuclear shadowing effects? 

It has been observed [33] that the quark distributions in a nucleus is depleted in the 
low region of x, and this is called nuclear shadowing. Our parton distributions [8] include 
nuclear shadowing corrections through the QCD evolution equations. To illustrate the 
shadowing effects in the multiplicity productions, we use the following scaled distributions 
by < Ncou >=< 0.5Npart >< z^ > to define the nuclear shadowing factor 
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The results in Au — Au collisions at ^/s = 200GeV (solid curve) and ^i = 5.5TeV 
(dashed curve) are presented in Fig.25. The calculations are stopped at larger rj, where 
the quark recombination mechanism and Fermi motion effects [30] become important. The 
77-dependence of the nuclear shadowing factor R presents a strong energy- dependence. We 
find that the shadowing corrections to the initial unintegrated gluon distribution cannot 
be negligible in any studies of the nuclear effects. 

(iv) Two component model or single component model? 

The early two component picture of hadron-hadron collisions assumes that a broad 
boost-invariant central plateau is separated by two energy-independent fragmentation 
regions [17]. However, the following experiments indicate that limiting fragmentation can 
be extended to beyond fragmentation ranges, and no evidence shows a boost-invariant 
central plateau [25]. It seems that a single-component mechanism dominates the hadron 
production. The KLN [6] and Albacete [7] are such single-component models. However, 
the dashed curves in Figs. 3 and 4, which are determined by the evolved unintegrated 
gluon distributions, are lower than the data, and it shows that the contributions from two 
components of production are necessary. 

(v) Saturation, it comes or not? 

Saturation is a limiting form of the shadowing modified gluon distribution. In this 
work, we do not consider the corrections of saturation, although shadowing is included. 
We compare our predictions with those of two saturation models. Kharzeeva, Levin, and 
Nardi [6] use a single production mechanism [i.e., the component I in Eqs. (12) and (14)] 
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to predict the heavy- ion colhsions at LHC. Instead of the two scale unintegrated gluon 
distribution in Eq. (1), KLN use a saturated integrated gluon distribution 

2x / a (Q^) ^Pti^ ~ ^)'^ Pt<Qs{x) 

to calculate the rapidity distributions. S is the inelastic cross section for the minimum 
bias multiplicity and /t is a normalization coefficient. A free parameter is dN/dr]\rj=o, 
which contains nonperturbative information. In the KLN model, it is fixed by the value 
of dNjiu_Au/dv\v=o ^^ V^ — 130 GeV and assumption 

< V{y/s) > c{y/s) = const. (22) 

Figure 26 compares our results in Fig. 17 with the KLN predictions [6]. 

Albacete [7] uses the same single production mechanism but with a different gluon dis- 
tribution, which is the solution of the BK equation in the form of single scale distribution. 
Albacete assumes 

< NpartiVs) >< T^iVs) > C{y/s) = COUSt., (23) 

and uses the value dNAu-Au/dri\r,=o at ^i = 200GeV to fixed this constant. Figure 27 
compares our results with the Albacete predictions [7]. 

Except for the contributions of quark recombination and different values of Trieff, the 
differences among three models originate from the different parton distributions, which 
obey different QCD evolution dynamics. Therefore, the observations of the multiplicity 
distributions are useful in identifying a true QCD evolution dynamics in the upcoming 
LHC data. 
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(vi) Validity of the factorization formula Eq. (1). We noted that the works in Ref. 
[34] presented a covariant gauge calculation, where the transverse momentum spectrum 
of the gluon is perturbatively determined by the final-state interactions of the gluon 
with the nucleons in the nucleus. In this case, the applications of the kt factorization 
formula (1) is unsatisfied. However, these works neglect the QCD evolution of gluons. 
Our aim is to test the predictions of the MD-DGLAP evolution equation in the LHC 
physics, where a physical gauge is used for the factorization of the evolution kernel. The 
above-mentioned final-state interactions are absent, and these effects are absorbed into 
the phenomeno logical fragmentation functions. Therefore, the application of Eq. (1) is 
reasonable in this work. 
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6 Summary 

We use two complementary production mechanisms: hard gluon-gluon fusion in the 
central rapidity region and soft quark recombination in the fragmentation region to study 
the particle multiplicity distributions in hadron-hadron collisions at high energies. We 
emphasize that a set of consistent integrated and unintegrated parton distributions, which 
are well defined in a broad kinematic range are necessary in such partonic model. For this 
reason, our parton distributions proposed in Ref. [8] are used. Based on the explanations 
of the present data, we predict the pseudorapidity densities of charged particles produced 
in p — p and central Pb — Pb collisions at LHC energies. We find that the limiting 
fragmentation hypothesis, which generally appears in the present data of hadron collisions 
is partly violated if the observations are across a wide range between the RHIC and 
LHC energies. An explanation about limiting fragmentation and its violation in partonic 
picture are given. We proposed that a possible QGP effect may deform the shape of 
the central plateau. In this work, we use an entirely phenomenological parametrization 
of existing data to fix the energy dependence of normalization factors of the hadronic 
distributions. This simple method may need modification after we obtain the LHC data 
at midrapidity. However, once this normalization is fixed, whole particle multiplicity 
distributions are completely determined by our parton distributions. The comparisons 
between our model and other partonic models are given. The differences in the predicted 
distributions in various models can help us to identify the true QCD dynamics in hadron 
collisions. 
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Fig. 1 Coefficient c{^/s) in Eq.(5) plotted as a function of energy ^/s using Eq.(13). 
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Fig. 2 Computed pseudo-rapidity distribution of charged particles in p — p(j)) collisions 
at various energies with rrieff = 770MeV and 6 = 0.7 (solid curves). The dashed curves 
take fUeff = 0. The data are taken from Ref. [20]. 
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Fig. 3 Solid curve in Fig. 2 at ^/s = 130GeV, where dashed and dotted curves are the 
contributions of gluon fusion mechanism I and quark recombination mechanism II. 
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Fig. 4 Similar to Fig. 3 but at ^s = 200GeV. 
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Fig. 5 Predicted pseudorapidity distribution of charged particles in p — p collisions at 
^/s = l.STeV. Dashed and dotted curves are the contributions of the gluon fusion model 
(I) and quark recombination model (II). 
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Fig. 6 Similar to Fig. 5 but at ^/s = 5.5TeV. 
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Fig. 7 Similar to Fig. 5 but at ^/s = 7TeV. 
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Fig. 8 Similar to Fig. 5 but at ^s = lOTeV. 
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Fig. 9 Similar to Fig. 5 but at ^/s = 14TeV. 
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Fig. 10 Comparisons of our predictions (solid curve) for pseudorapidity distribution of 
charged particles in p — p collisions at ^/s = 14TeV with that of KLN predictions (dashed 
curve) [6]. 
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Fig. 11 Comparisons of our predictions (solid curve) for pseudo-rapidity distribution 
of charged particles in p — p collisions at ^/s = 14TeV with that of PYTHIA [22] and 
PHOJET models [23]. 
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Fig. 12 Comparisons of our predictions (solid curve) for pseudorapidity distribution of 
charged particles in p — p collisions at ^/s = 14TeV with that of the UrQMD model 
(dashed curve) [24]. 
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Fig. 13 Pseudo-rapidity density of charged particles produced in Au — Au collisions with 
0—6% central cut at ^/s = 130 GeV. Data are taken from [29]. The solid curve corresponds 
to TTLeff = 770MeV. Dashed and dotted curves are the contributions of the gluon-gluon 
fusion and quark recombination, respectively. 
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Fig. 14 Same as Fig. 13 but for y^ = 200GeV. 
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Fig. 15 Pseudorapidity density of charged particles produced in Au — Au collisions with 
— 6% central cut at ^/s = ISOGeV with different parameters irieff = 770MeV (solid 
curve), 400MeV (dashed curve), and (dotted curve). Data with errors are taken from 
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Fig. 16 Same as Fig. 15 but at y/s = 200GeV. 
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Fig. 17 Predictions of pseudo-rapidity density of charged particles produced in Ph — Pb 
collisions with — 6% central cut at ^/s = 5.5 TeV. Dashed and dotted curves are the 
contributions of the gluon-gluon fusion and quark recombination, respectively 
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Fig. 18 Shifted and scaled pseudorapidity distribution of charged particles produced in 
Au — Au collisions with — 6% central cut at ^/s = 130GeV (dashed curve) and 200GeV 
(solid curve). Results exhibit limiting fragmentation at rj' > —1.5. 
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Fig. 19 Similar to Fig. 18 but with contributions from gluon-gluon fusion mechanism. The 



same hmiting fragmentation appears. 
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Fig. 20 Comparisons of shifted and scaled pseudorapidity distribution of charged particles 
produced in Pb — Pb collisions with — 6% central cut at ^/s = 5.5TeV (dotted curve) 
with the curves of Fig. 18. Results exhibit deviations from limiting fragmentation at 

7]' <0 
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Fig. 21 Similar to Fig.20 but with contributions from the gluon-gluon fusion mechanism. 
Similar violations of limiting fragmentation appear. 
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Fig. 22 Predicted dNj^^N/drj' at 0—6% central vs the shifted pseudo-rapidity rj' = i]—yi,eam 
in a range of energies corresponding to Figs. 3-9 {^/s = 130GeV-14TeV from bottom to 
top). The results exhibit longitudinal scaling and a small violation as in Fig. 21. 
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Fig. 23 Kinematical ranges of Xi/2 with different energies ^/s and pt = 0.5GeV in Eq.(2): 
Solid curves: y/s = 5.5TeV; Dashed-dotted curves: ^i = 200GeV and dotted curves: 
y/s = 130GeV. Vertical dashed curves correspond to (from left to right) y' = at 
v/s = ISOGeV, 200GeV and 5.5GeV. 
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Fig. 24 Deformation of central rapidity plateau with decreasing parameter rrieff in Pb—Pb 
collisions at y/s = 5.5 TeV. Solid, dashed, and dotted curves correspond to m^ff = 
400MeV, 200MeV and 0. 
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Fig. 25 Nuclear shadowing factor R in Eq.(20) at different energies ^/s, where the pseu- 
dorapidity distributions in Au — Au colhsions are scaled by Ncou- 
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Fig. 26 Comparison of our predictions in Fig. 17 with the resuhs of the KLN model [6] 
(dotted curve), where mefj = 500 MeV in the KLN modeL 
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Fig. 27 Comparison of our predictions in Fig. 17 with the results of the Albacete model 
[7] (dotted curve), where meff = 250MeV in the Albacete model. 
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